Forced MHD turbulence in three dimensions using Taylor-Green symmetries 
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We examine the scaling laws of MHD turbulence for three different types of forcing functions 
and imposing at all times the four-fold symmetries of the Taylor-Green (TG) vortex generalized to 
MHD; no uniform magnetic field is present and the magnetic Prandtl number is equal to unity. We 
also include a forcing in the induction equation, and we take the three configurations studied in the 
decaying case in [E. Lee et al. Phys. Rev.E 81, 016318 (2010)]. To that effect, we employ direct 
numerical simulations up to an equivalent resolution of 2048 3 grid points. We find that, similarly 
to the case when the forcing is absent, different spectral indices for the total energy spectrum 
emerge, corresponding to either a Kolmogorov law, an Iroshnikov-Kraichnan law that arises from 
the interactions of turbulent eddies and Alfven waves, or to weak turbulence when the large-scale 
magnetic field is strong. We also examine the inertial range dynamics in terms of the ratios of 
kinetic to magnetic energy, and of the turn-over time to the Alfven time, and analyze the temporal 
variations of these quasi-equilibria. 



PACS numbers: 47.27.-i, 47.65.-d 

I. INTRODUCTION 

Turbulence is a common feature of a variety of flows, 
from engineering to geophysics and astrophysics. It re- 
mains unsolved, due in part to a lack of statistical theory 
on how to deal with a very large number of modes in- 
teracting nonlinearly, and competing with waves. At the 
moderate Reynolds numbers that are achievable today 
numerically in three space dimensions on uniform grids 
of at most 4096 3 points, one follows accurately the tem- 
poral evolution of in excess of 64 billion modes, leading to 
the creation of myriads of vortex filaments. When cou- 
pling to a magnetic field, in the magnetohydrodynamic 
(MHD) approximation regime for velocities small com- 
pared to the speed of light, one observes current and 
vorticity sheets, that are found to roll-up for sufficiently 
high Reynolds numbers [T] . 

One question concerns the universality or not of the 
scaling laws of turbulent flows. There has been much 
debate concerning this point, in particular in the MHD 
community as well as when dealing with the dynamics of 
the atmosphere and the oceans: one way to phrase the 
question is to ask wether the presence of waves will affect 
the energy distribution among modes, inertial waves in 
the rotating case with solid body rotation, gravity waves 
for stratified turbulence, Alfven waves in MHD, acoustic 
waves when the condition of incompressibility is removed, 
as is necessary in the interstellar medium where super- 
sonic flows are routinely observed. The answer is unam- 
biguous in the regime of weak turbulence when the ratio 



of characteristic times (wave period and eddy turnover 
time) is small; this small parameter allows for a natu- 
ral closure to the statistical problem and constant-flux 
(as well as zero flux) solutions can be found in terms of 
power laws as a function of anisotropic wave numbers, 
the anisotropy arising from the imposition of an exter- 
nal agent (uniform rotation, gravity, or magnetic field), 
and to the anisotropic dispersion relations [2J. But this 
weak-turbulence regime is non-uniform in scale, simply 
because the variation with scale of the wave period tw 
and of the eddy turn-over time tnl are different; hence, 
there exists a scale at which these two timescales are 
equal and the weak turbulence regime breaks down. For 
stratified flows, this is called the Ozmidov length scale, 
and for the rotating case, the Zeman scale. Note that for 
MHD, the situation is different: for stratified and rotat- 
ing flows, at scales smaller than the Ozmidov or Zeman 
scales, isotropy and a classical Kolmogorv scaling is likely 
to recover (3], whereas it does not in MHD. In fact, one 
could argue the opposite: isotropy can prevail at large 
scales where the effect of the large-scale magnetic field is 
purely local and its amplitude is comparable to that of 
the modes it is interacting with, whereas 

In MHD, another hypothesis has been put forward to 
understand the dynamical exchanges in a phenomcnologi- 
cal way, that of an equality between the two characteristic 
time scales, an equality that would hold throughout the 
inertial range [14 . This hypothesis leads to a Kolmogorov 
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spectrum E(k±) ~ k x (hereafter Kp41), expressed in 
terms of k±, where the direction refers to that of the 
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external agent, here a uniform magnetic field. It was 
found in [5] that the same hypothesis can also lead to an 
Iroshnikov-Kraichnan spectrum (IK hereafter) or a weak 
turbulence spectrum, E(kx,k») ~ kj 2 f(k^) (WT here- 
after), on the simple basis that Tw(k)/rNL(k) = r T (k) 
can be constant but not necessarily equal to unity: the 
different regimes appear in that light as emerging from 
a different rate at which energy is exchanged between 
its kinetic and magnetic modes. All these spectra have 
been observed in direct numerical simulations (DNS) of 
three-dimensional (3D) MHD turbulence. In one par- 
ticular case, identical velocity fields are used as ini- 
tial conditions, with comparable invariants (total en- 
ergy E tot = 0.25 with E v (t = 0) = E b (t = 0) where 
E v (, are the kinetic and magnetic energy respectively, 
with strictly zero magnetic helicity and negligible cross- 
correlation between the velocity and the magnetic field) 
[6j. It is the purpose of this paper to pursue the work 
done in [BJ, extending it to the statistically steady case 
in the presence of forcing. 

In the next section, we write equations, initial con- 
ditions and forcing; the results are given for a high- 



resolution run in ( III and in £ IV we compare the evo- 



lution of three magnetic configurations, at resolutions of 
1024 3 grid points. Finally, qVjis the conclusion. 



II. THE EQUATIONS AND THE NUMERICAL 
SET-UP 



The MHD equations for an incompressible fluid with 
v and b respectively the velocity and magnetic fields in 
Alfvenic units are: 

+ v • Vv = VP + j x b + iA7 2 v + F v , (1) 

at po 



db 

~dt 



V x (v x b) + ?yV 2 b + F M ; 



(2) 



po = 1 is the (uniform) density, and b is dimensionally a 
velocity as well, the Alfven velocity; V is the total pres- 
sure, V ■ v = V • b = 0, and v and r\ are respectively 
the kinematic viscosity and magnetic diffusivity; we take 
v = r\ (unit magnetic Prandtl number). Finally, Fv,Fm 
are forcing terms introduced both in the momentum and 
in the induction equations. In principle, above a criti- 
cal Reynolds number R M , a dynamo mechanism sets in 
whereby sufficient magnetic excitation is produced at all 
scales. For the Taylor-Green flow defined below it was 
shown in reference [7] that R M depends very strongly on 
the imposed symmetries. In addition, when imposing all 
symmetries at all times, it was shown in reference [8] that 
R M is very high (of the order of 1000). In this work, like 
in reference [BJ , we focus of the fully symmetric problem 
in order to maximize the available Reynolds number, and 
thus the maximum resolution (see discussion following 
Eq. (14) below). Hence, we chose to force the induction 



that, by simply breaking the symmetry of the initial con- 
ditions and using a general code, this critical parameter 
is lowered by more than one order of magnitude, but with 
a substantially costlier computation, by a factor of 32 [5] . 

The energy E to t, the cross helicity Hq and the mag- 
netic helicity Hm are defined as 



-Etot = 



■) + E b 
(v-b) 



H 



M 



b 2 )/2 
= (A • b) 



(3) 
(4) 



where b = V x A and A is the magnetic potential. In 
the ideal case (v = rj = 0) and without forcing (Fy — 
Fm = 0) note that these quantities (E tot ,Hc and Hm) 
are all conserved. Relative helicities measure the relative 
alignment of vectors, independently of their amplitudes, 
pv = [cosv,o;] , pc — cos [v, b] , pm — cos [A, b], with 
lu = V x u the vorticity (the kinetic helicity Hy = (u-oj) 
is an invariant when b = 0). 

Considering a flow which is 27r-periodic in all spatial 
dimensions, the kinematic Reynolds number Re and the 
magnetic Reynolds number R D are defined as 



Re = 



Lv ri 



Rh — 



(5) 



where the root-mean square velocity is t> rms = y^2E v /3 
and the characteristic length L is defined by 



L = 2tt 



J2 k E v (k,t)dk ' 



(6) 



where the kinetic energy spectrum E v (k,t) (such 
that E v (t) = *^2 k E v (k,t)) is obtained by summing 
i|u(k',i)| 2 on the spherical shells A;— 1/2 < |k'| < fc+1/2 
(u(k) is the Fourier transform of the velocity). Anal- 
ogously, the magnetic energy spectrum is denoted by 
E b (k,t) and verifies E b (t) = J2k E b {k,t). 

We now turn to the definition of the external driving 
volumic forces Fy,M in Q-fl^]) which balance the total 
energy dissipation and allow to reach a statistically sta- 
tionary state. Following reference 0, we force the system 
by setting in ([!]) 

Tp f _,TG 



(7) 



where v TG is the Taylor- Green vortex [9] given by 



,TG 



(sin(x) cos(y) cos(z), — cos(a;) sin(y) cos(,z), 0), 



(8) 

and f v is always set to the value 1/16. The force Fm is 
determined in a similar way but using instead of the TG 
velocity mode, the three magnetic field modes studied 
in the decaying MHD runs of reference [6j. These three 
modes were: 



cos x sin y sin z 
sin x cos y sin z 
2 sin x sin y cos z . 



b 1 = 6q I sin x cos y sin z 



cos 2x sin 2y sin 2z 
by — sin 2x cos 2y sin 2z 



(9) 



(10) 



equation as a way to mimic the dynamo itself. We note 
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sin 2x cos 2y cos 2z 
cos 2x sin 2y cos 2z 
-2 cos 2x cos 2y sin 2z , 



(11) 



The labels /, A and C stand respectively for insulat- 
ing, alternate insulating and conducting boundary con- 
ditions for the current when considering its orientation 
with respect to the wall of the so-called fundamental box 
(see [Ifl]). The coefficients b^, b$ and b§ are such that, 
for all cases, = 1/16. Thus, all computations have 
equal initial kinetic and magnetic energy at t=0, with 
■Etot = E v + Ef, = 1/8. 

Like in Eq. ^ , the amplitudes of the forcing are cho- 
sen as 



Fm — fbb l ' A ' C , 



(12) 



and the values of /& are given below in Table [TJ 

We shall also examine the behavior of the spectral ra- 
tios of time-scales R T and of modal energies Re defined 
respectively as 

Rr(k) = T NL (k)/T A (k), R E (k) = E b {k)/E v (k) , (13) 

with TNh{k) = [ku(k)]^ 1 and r^fc) = [fcS ] _1 , -Bo being 
defined here as the magnetic field in the gravest mode. 

Because of the symmetries of the TG vortex extended 
to MHD, all fields can be represented in Fourier space as: 

oo 

v x( r )t) = u x (to, n,p,t) sin mx cos ny cos pz, ( 14) 

m,n,p— 
oo 

v y (r,t) = u y (m, n,p, t) cosmx sinny cospz,(15) 

m,n,p— 
oo 

v z (r,t) = U z (m, n,p, t) cos mx cosny sinpz,(16) 

m,n,p— 

where u(to, n,p, t) is equal to zero unless the three in- 
tegers m,n,p are either all even or all odd. Thus, in 
spectral space, the following symmetry is fulfilled: 

u x (m,n,p,t) = (-l) r u y (n,m,p,t) , (17) 
u z (m,n,p,t) = (-l) r u z (n,m,p,i) , (18) 

with r = 1 if m, n,p are all odd, and r — 2 if m,n,p are 
all even. 



Relations (17) and (18) allow one to only compute v x 
and v z . Moreover, the decomposition (14) on cither even 



or odd integers leads to a gain of a factor of 32 in memory 
and CPU time compared with the general case of Fourier 
transforms with the same scale separation, or k max /k m i n , 
with respectively k m i n = 1 for a box of length 2ir and 
kmax = N/3 with N the number of grid points per di- 
mension, using a standard 2/3 de-aliasing rule. The code 
is pseudo-spectral, with a fourth-order Runge-Kutta tem- 
poral scheme and with periodic boundary conditions. All 
previous symmetry relations are implemented to speed 
up calculations. The code is parallelized up to ~ 98,000 
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TABLE I: List of runs and parameters. Resolution JV, final 
time Tf and time step dt, viscosity and magnetic diffusivity 
v — r), forcing parameter /j, and Reynolds number Re (see 
Eqs. flppl andp])). 



processors on grids of up to 8196 3 points, using a hybrid 
(MPI-Opcn-MP) algorithm which becomes advantageous 
at high resolution |llj . Grids used in this work have the 
equivalent resolution of 1024 3 and 2048 3 points. Run 
parameters are summarized in Table [I] 



III. RESULTS FOR THE C RUNS 

The first run on which we report is the one at the 
highest Reynolds number (and the highest resolution). 
We give in Fig. [T]a, top, the temporal evolution of the 
kinetic, magnetic and total energy. The thick lines are 




FIG. 1: (Color online) a: Temporal evolution of kinetic 
(blue), magnetic (red) and total (black) energy for Runs CI 
and C2 (see Table |B; b: Magnetic to kinetic energy ratio as 
a function of time for the same runqs; this ratio seems to 
saturate at a value ~ 0.55. In both figures, longer temporal 
evolutions up to t ~ 65, are given for the run on a grid of 
1024 3 points. 

for Run CI and the dashed lines for run C2 at lower 
resolution (see Table [I]) . The ratio of magnetic to kinetic 
energy Re is given in Fig. [l]b. The total energy saturates 
for times larger than ~ 20, at ~ 0.12, and the global 
energetic ratio is between Re ~ 0.3 and Re ~ 0.6, with 
a tendency toward growth. 

We performed a temporal average of the total energy 
spectrum from t = 17 A to t = 22.5 for runs CI (red, 
solid line) and C2 (blue, dashed line) which we present 
in Fig. |2|a, compensated by fc +3 / 2 . As in the decay 
case presented in [6], the best fit is for E to t(k) ~ fc~ 3 / 2 , 
i.e. an Iroshnikov-Kraichnan law. We give further below 
the temporal evolution of the instantaneous value of the 
spectral index in Fig{7]b (bottom), when comparing it 
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FIG. 2: (Color online) For runs CI and C2 (see Table [J, a) 
fc 3//2 -compensated total energy spectrum temporally averaged 
in the interval t = 17 — 22.5 for run CI and t = 20 — 68 for run 
C2. b) Spectral ratios of energy -Re(&) (blue) and time-scales 
R T {k) (red), as defined in Eq. (13 1. Note the constancy of 
the former in the inertial range. 



for several runs (see next Section). 

In Fig. [2jb is shown the ratio of the spectra of the 
turn-over time to the Alfven time, R Tl defined in Eq. 



(13) and averaged in the same temporal interval as the 
energy spectrum. As can be seen from Fig|2]b, there is 
a systematic increase of this ratio in the inertial range, 
contrary to the hypothesis of critical balance advocated 
in [4] (see also Q2])- Rather, it is the magnetic to 
kinetic spectra ratio which remains remarkably constant 
throughout the inertial range, as displayed in Fig[2jb, 
with as often, a slight excess of magnetic energy, except 
at the gravest mode that dominates the global energetics; 
this confirms the earlier findings of the decay case [5] , as 
well as those in numerous other numerical simulations 
(see e.g. [13]). 

When examining the behavior of the C-iiow at lower 
resolution on a grid of 1024 3 points, we observe that the 
results are in agreement with these conclusions; the lower 
resolution simply allows us to compute for longer times, 
leading to a better temporal averaging. Nevertheless, 
there may be a trend toward the energetic ratio to in- 
crease at later times (see Fig. [I]). 



IV. COMPARATIVE RESULTS FOR THE 
THREE FORCING FUNCTIONS 

The most striking result of the computations per- 
formed in the decaying case presented in [S] is that dif- 
ferent initial conditions for the magnetic field only, but 
with the same global invariants, led to different spectral 
inertial indices. Will the same occur in the presence of 
forcing? This is what we are now investigating. We thus 
address now the question of the scaling of the two other 
configurations studied in [6, in the decaying case, namely 
the so-called A- and I- magnetic configurations. 



A. The A run 

The forcing with the A configuration leads the system 
to reach a statistically stationary state, both globally for 
Etot (shown in Fig. [3ja) and in its kinetic to magnetic 
energy ratio (as displayed in Fig. |3]b), with a value of 
that ratio slightly above unity, as often observed in the 
Solar Wind. The total energy spectrum for this run is 
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FIG. 3: (Color online) a)Temporal evolution of kinetic (blue), 
magnetic (red) and total (black) energy for Run A (see Table 
llj. b) Temporal evolution of the ratio Et/E v for the same 
run. Equipartition is almost reached in that run, with a slight 
excess of magnetic energy. 

rather steep, with the best fit corresponding to a fc~ 2 
law (see Fig. |4]a), i.e. a law corresponding to weak 
turbulence. This is in contrast to [6]: in the decaying 
case, this flow had a spectral index close to —5/3. Note 
that the wave-turbulence behavior in the present case is 
consistent with the ratio of energy spectra and the ratio 
of time scales displayed on Fig.|I]b. 
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FIG. 4: (Color online) Run A (see Table [Tj; all spectra are 
averaged in time in the interval t = 25 — 52 and computed 
by summing on spectral shells of width Ak = 3 to reduce 
oscillations, a) fc 2 -compensated total energy spectrum, b) 
Ratios of energy spectra (red), again constant in the inertial 
range, and ratio of time-scales (blue). 



In Fig. [7]b below the temporal evolution of the spec- 
tral indices of all the runs are compared. Note that the 
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spectral index of the total energy spectrum of the A-flow 
varies substantially over time. At this point, it is diffi- 
cult to decide which power law is best followed. Initial 
behavior seems consistent with the unforced —5/3 value 
although the inertial law is steepest at later times, and 
is thus more in favor of a weak turbulence spectrum. 
However note that steep structures, such as sharp and 
isolated current and vorticity sheets, can also lead to a 
"shock" like spectrum. In this context, see the visualiza- 
tions presented below in Fig. [8] at the end of the present 
section. 



B. The I runs 

Let us now examine the dynamics of the magnetic I 
configuration (see definition in ([9])). The temporal evo- 
lutions of the energies is displayed in Fig(5]a. The total 















E « E b Etot 



















5 


10 


15 20 


25 

t 


30 


35 40 45 50 




























| E b /E v | 

















°0 5 10 15 20 25 30 35 40 45 50 



t 




FIG. 6: (Color online) Run I (see Table|T]|; all spectra are av- 
eraged in time in the interval t = 30—52 , a) fc 5 / 3 -compensated 
total energy spectrum, b) Ratios of energy spectra (blue) and 
time-scales (red). 



As already observed in [B] in the decaying case, it may 
be that spectral indices vary with time, as the ratio of 
kinetic to magnetic energy varies as well. Also, it is quite 
difficult to distinguish between spectral indices that are 
quite close and this flow seems more undecided than for 
the other two flows we study in this paper. 

To summarize, in FigjTJa we present the (uncompen- 
sated) spectra for the three runs (averaged over time). 
We also show in Fig. [Tja (bottom) the temporal variation 
of spectral indices for the three runs. 



FIG. 5: (Color online) a)Temporal evolution of kinetic 
(blue), magnetic (red )and total (black) energy for Run / 
(see Table III), b) Temporal evolution of the ratio Et/E v for 
the same run. 



energy seems to reach a quasi steady state although an 
increase of the ratio of magnetic to kinetic energy is ob- 
served in Fig(5]b. Thus, a stationary state is not reached 
for that ratio, with marked global oscillations. 

The total energy spectrum (compensated by fc 5//3 ) and 
the time-scales ratio and energy spectra ratio are shown 
in [6] A correlation between the oscillatory growth of the 
ratio of magnetic to kinetic energy displayed in Figj5]b 
and the spectral index of run I displayed as the middle 
curve in Fig(7]b is apparent. Indeed upward trends in the 
spectral index, moving towards weak turbulence, corre- 
spond to growth of magnetic energy at the expense of 
the kinetic energy. Observe in Fig|6]b that again we ob- 
tain close to equipartition in the medium range of scales, 
with another increase of magnetic energy in the smallest 
scales. The variation of time-scales is dominated by the 
growth in the gravest mode of the magnetic energy since 
the Alfven time is controlled by the large scale magnetic 
field Bq defined as equivalent to E},(k = 2). 



a) 10" 

10" 
W 10" 







k- 5 ' 3 


-^Run C1 
-^Run A 
-»-Run I 

















10 k 10 10 

RUN C1 RUN C2 RUN A RUN It 




FIG. 7: (Color online) a) Total energy spectrum for the 
runs C2, A and /. Black dashed lines display the different 
fits, b) Temporal evolution of spectral indices for all runs 
(see inset). Note that the differentiation between the three 
magnetic configurations is clear at all times, even if the exact 
value of the indices vary with time. 



For completeness and in order to illustrate the 
physical-space distribution of magnetic energies in the 
the respective runs, we have performed some visualiza- 




tions using the VAPOR software [14]. It is apparent on 
FigjH] that large-scale coherent structures are present in 
the I and A flow. In contrast, the C-flow seems more 
isotropic. As discussed above it is a possibility that steep 
structures, such as sharp and isolated current and vor- 
ticity sheets, can also lead to scaling in energy spectrum. 
This certainly looks as a possibility in the cases of I and 
A flows. Further studies of these structures are left for 
future work. 



V. CONCLUSIONS 

In this paper, we have extended the analysis of non- 
universality of MHD spectra from the decaying case per- 
formed in |6] to the forced case. We confirm the previ- 
ous results, with either IK, K41 or WT spectra emerging 
on average, when using Taylor-Green forcing, including 
in the induction equation, although temporal variations 
may be occurring. Note that lack of universality in MHD 
has already been found by other authors, in the context 
of heating the solar corona [T5lU7| . or in the presence 
of strong correlations between the velocity and the mag- 
netic field [i"8l420j . We also confirm that these different 
scalings are linked to the magnetic energy content in the 
gravest mode and that, at least for these flows in which 
the four-fold symmetries of the TG vortex are imposed 
at all time, quasi-equipartition between the kinetic and 
magnetic energy obtains in the inertial range, with a vari- 
ation of the ratio of the turn-over time to the Alfven time 
consistent with the inertial index of the energy spectrum, 
as in the decaying case and contrary to the hypothesis 
made in [J]. 

It should be noted that, even though there is no im- 
posed magnetic field in these computations and thus no 
imposed anisotropy, an extension of this work could be 



to analyze the data in terms of anisotropic scaling with 
respect to a locally-defined quasi-uniform field, averag- 
ing the induction in a sphere of diameter the integral 
scale, as done for example in [5T] (see and references 
therein for anisotropic scaling in MHD). It is also well 
known that, in the atmosphere and the oceans, different 
spectra may emerge according to the relative strength 
of the stratification, the rotation and the forcing, as for 
example the Garrett-Munk versus Phillips spectra [23] . 

There may be periods of evolution when the flow tends 
to one regime, and at other times to another regime. It 
was already observed in [6] that for late times, the dis- 
tinction between the K41 and IK regimes became diffi- 
cult to make but of course the Reynolds number by then 
had decreased substantially. It is known that there are 
long-time fluctuations in most turbulent flows ([24] and 
references therein) and this could lead to alternate ex- 
changes of behavior, as already observed in |25j in two- 
dimensional MHD, with turbulent periods of the order of 
one hundred turn-over times. These long-time fluctua- 
tions could lead to long-time fluctuations in spectral in- 
dices as well. Such behavior can be attributed to the last- 
ing effects of nonlocal interactions between widely sepa- 
rated scales [3S] H7J as observed in high-resolutions DNS 
of MHD turbulence. 

It is not known whether these results will stand out 
in the limit of infinite Reynolds numbers, and higher 
Reynolds number computations will have to be per- 
formed in order to confirm the results presented in this 
paper. There are other venues that can be taken as well: 
it is well-known that magnetic helicity and cross-helicity 
play essential roles in the dynamics oh MHD turbulence, 
and yet they arc quenched by the symmetries in the 
present approach where symmetries are enforced at all 
times. It is already documented that, for the case of 
long time dynamics and perturbing the three initial con- 
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ditions studied in this paper, vastly different regimes can 
be reached, with in particular the ratio of kinetic to mag- 
netic energy varying in a large range |13j . This behavior 
can be understood in terms of minimization of energy 
subjected to the constraints of the invariance of He and 
Hm [2H1 ES]. This points out to the possibly essential 
role played by the imposed four-fold symmetries. Such 
symmetries could be broken in part, as performed in [7], 
leading to different modes growing in the dynamo regime. 
A similar approach could be taken for the present prob- 
lem of lack of universality, still allowing for some sav- 
ings in computer resources compared with full-fledged 
MHD computations which might otherwise have to be 
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